####### Robustness tests ##############

# This script generates tables and figures showing robustness tests.
# All tables are displayed in Appendix B.
# Figure 6a and 6b were displayed in the main text.

# This script produces:
# Figures 6a, 6b, B1
# Tables B1 to B35
# Figures B1 and B2


# Session info
# R version 4.1.0 (2021-05-18)
# Platform: x86_64-apple-darwin17.0 (64-bit)
# Running under: macOS 13.1


# Figure 5a ---------------------------------------------------------------

threshold_list <- c(8,15,22,30)
for (threshold in threshold_list) {
  
  load(paste0("data/navco_week_replication",threshold,".Rda"))
  
  
  nb4 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  #time
                  diff + diff2  +protests_past + pr_nr_c + 
                  factor(gid_year),
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))  
  
  assign(paste0("nb4_", threshold, sep=""), nb4)
  
}

### Coefficient plot
m1_df <- tidy(nb4_8) %>% 
  filter(term == "con_sum") %>% 
  mutate(model = "7 days")

m2_df <- tidy(nb4_15) %>% 
  filter(term == "con_sum") %>% 
  mutate(model = "14 days")

m3_df <- tidy(nb4_22) %>% 
  filter(term == "con_sum") %>% 
  mutate(model = "21 days")
m4_df <- tidy(nb4_30) %>% 
  filter(term == "con_sum") %>% 
  mutate(model = "30 days")


many_models1 <- rbind(m4_df, m3_df, m2_df, m1_df) %>% 
  mutate(signicant = as.character(ifelse(p.value < 0.05, 1, 0)))

dwplot(many_models1, 
       vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
       dot_args = list(aes(shape = signicant))) %>% 
  relabel_predictors(c(con_sum = "Concession (sum)",                       
                       size_log = "Number of protesters (log)",
                       protesterviolence = "Protester violence",
                       goal_policy = "Demand: policy change",
                       goal_reform = "Demand: institutional reform",
                       goal_regime = "Demand: regime change",
                       state_violence = "State violence")) +
  scale_shape_discrete(name = "",
                       labels = c("p-value < 0.05")) +
  theme_lh + theme(legend.box = "vertical") +
  scale_colour_grey(start = .7, end = 0, name = "Time window") +
  xlab("Coefficient estimate") + 
  xlim(-.05, .5) +
  ylab("")
ggsave("output/figures/Figure5a.pdf", scale = 1, dpi = "retina", width = 19, height = 10, units = "cm")

# Figure 5b:  --------------------------------------------------------
threshold_list <- c(8,15,22,30)
for (threshold in threshold_list) {
  threshold1 = threshold -1
  load(paste0("data/navco_week_national_replication",threshold,".Rda"))
  
  nb4 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime + state_violence +  
                  #time
                  diff + diff2 + protests_past +
                  factor(country_year),
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))

assign(paste0("nb4_", threshold, sep=""), nb4)
}

m1_df <- tidy(nb4_8) %>% 
  filter(term == "con_sum") %>% 
  mutate(model = "7 days")

m2_df <- tidy(nb4_15) %>% 
  filter(term == "con_sum") %>% 
  mutate(model = "14 days")

m3_df <- tidy(nb4_22) %>% 
  filter(term == "con_sum") %>% 
  mutate(model = "21 days")
m4_df <- tidy(nb4_30) %>% 
  filter(term == "con_sum") %>% 
  mutate(model = "30 days")


many_models2 <- rbind(m4_df, m3_df, m2_df, m1_df) %>% 
  mutate(signicant = as.character(ifelse(p.value > 0.05, 1, 0)))

dwplot(many_models2, 
       vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
       dot_args = list(aes(shape = signicant))) %>% 
  relabel_predictors(c(con_sum = "Concession (sum)",                       
                       size_log = "Number of protesters (log)",
                       protesterviolence = "Protester violence",
                       goal_policy = "Demand: policy change",
                       goal_reform = "Demand: institutional reform",
                       goal_regime = "Demand: regime change",
                       state_violence = "State violence")) +
  scale_shape_discrete(name = "",
                       labels = c("p-value < 0.05", "p-value > 0.1")) +
  theme_lh + theme(legend.box = "vertical") +
  scale_colour_grey(start = .7, end = 0, name = "Time window") +
  xlab("Coefficient estimate") + 
  xlim(-.05, .5) +
  ylab("")
ggsave("output/figures/Figure5b.pdf", scale = 1, dpi = "retina", width = 19, height = 10, units = "cm")


# Figure B1 ----------------------------------------------------------------
# Event study plot  
# add variables for Difference-in-differences analysis
threshold = 8
load(paste0("data/navco_week_replication",threshold,".Rda"))
nav_did <- navco_week %>% 
  group_by(move_id) %>% 
  mutate(con_first = ifelse(row_number() == 1, con_dummy, cumsum(con_dummy) + lag(con_dummy)),
         con_first = ifelse(con_first == 1, 1, 0),
         t_week = ifelse(sum(con_dummy) > 0, week[con_first == 1], 10000),
         treat = ifelse(sum(con_first) > 0, 1, 0)) %>% 
  ungroup() %>% 
  mutate(time_to_treat =  week - t_week,
         time_to_treat = ifelse(time_to_treat < -1000, 0, time_to_treat),
         post = ifelse(time_to_treat > -1, 1, 0),
         post_t = post*treat)

mod_twfe = feols(protests_future~ i(time_to_treat, treat, ref = -1)  +
                   # protest characteristic
                   size_log + protesterviolence + goal_policy + goal_reform + goal_regime + state_violence +
                   # time
                   pr_nr_c + diff + diff2 + protests_past|  
                   gid + week_calendar,                             
                 cluster = ~gid,                          
                 data = nav_did)

pdf("output/figures/FigureB1.pdf", width = 5, height = 4)
iplot(mod_twfe, xlim = c(-3,9), xlab = "Weeks before/after 1st concession", main = "Grid-level")
dev.off()

# Tables B1 to B3: other thresholds --------------------------------------------------------
threshold_list <- c(15,22,30)
for (threshold in threshold_list) {
  threshold1 = threshold -1
  
  load(paste0("data/navco_week_replication",threshold,".Rda"))
  
  nb1 <- glm.nb(protests_future ~ con_sum,
                data=navco_week)
  
  nb2 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb3 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  v2x_clphy + v2x_polyarchy +couprisk + nlights_calib_mean + gdp_log + pop_log +
                  #time
                  diff + diff2 + protests_past +pr_nr_c +
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb4 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence +
                  #time
                  diff + diff2 + protests_past +pr_nr_c +
                  factor(gid_year),
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))  
  
  # Table
  clustered_se <- list(
    sqrt(diag(vcovCL(nb1, vcov = vcovCL(nb1, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb2, vcov = vcovCL(nb2, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb3, vcov = vcovCL(nb3, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb4, vcov = vcovCL(nb4, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T)))
  )
  
  
  stargazer(nb1, nb2,nb3,nb4,
            se = clustered_se,
            type = "latex",
            float = F,
            float.env = "table",
            no.space = T,
            dep.var.caption =paste0("Sum of protests in the next ", threshold1, " days"),
            dep.var.labels = c(""),
            model.numbers = F,
            column.labels = c("1", "2", "3", "4"),
            title = paste0("Concessions and protests in the next ", threshold1, " days. Negative binomial models"), # is not written in latex if float=F
            notes.append = T,
            notes = c("Standard errors are clustered at the level of grid cells."),
            covariate.labels = c("Concession (sum)",
                                 "Number of protesters (sum, log)",
                                 "Protester violence (sum)",
                                 "Demand: Policy change",
                                 "Demand: Institutional reforms",
                                 "Demand: Regime change",
                                 "State violence (sum)",
                                 "Physical violence index",
                                 "Democracy index",
                                 "Coup threat",
                                 "Nightlights",
                                 "GDP pc (lag and log)",
                                 "Population (lag and log)",
                                 "Days since last protest",
                                 "Days since last protest (squared)",
                                 "Protest (lag)",
                                 "Protests in the rest of the country"),
            omit= c('year', 'gid', 'country'),
            omit.stat = c("logrank", "wald", "ll", "theta"),
            add.lines = list(c("Grid dummies", "No", "Yes","Yes", "Yes"),
                             c("Year dummies", "No", "Yes","Yes", "Yes"),
                             c("Grid-year interaction", "No", "No", "No", "Yes"),
                             c("Number of grids", "169", "169", "169", "169")),
            label = "table:main-negbin-national",
            out = paste0("output/tables/Appendix_Table_varyingTime",threshold,".tex"))
}


# Tables B4 to B7: country level ----------------------------------------------------------
threshold_list <- c(8,15,22,30)
for (threshold in threshold_list) {
  threshold1 = threshold -1
  load(paste0("data/navco_week_national_replication",threshold,".Rda"))
  
  nb1 <- glm.nb(protests_future ~ con_sum,
                data=navco_week)
  
  nb2 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  country_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb3 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  v2x_clphy + v2x_polyarchy +couprisk+ gdp_log + pop_log +
                  #time
                  diff + diff2 + protests_past  +
                  country_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb4 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence +  
                  #time
                  diff + diff2 + protests_past +
                  factor(country_year),
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  # Table
  clustered_se <- list(
    sqrt(diag(vcovCL(nb1, vcov = vcovCL(nb1, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb2, vcov = vcovCL(nb2, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb3, vcov = vcovCL(nb3, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb4, vcov = vcovCL(nb4, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T)))
  )
  
  
  stargazer(nb1, nb2,nb3,nb4,
            se = clustered_se,
            type = "latex",
            float = F,
            float.env = "table",
            no.space = T,
            dep.var.caption =paste0("Sum of protests in the next ", threshold1, " days"),
            dep.var.labels = c(""),
            model.numbers = F,
            column.labels = c("1", "2", "3", "4"),
            title = paste0("Concessions and protests in the next ", threshold1, " days. Negative binomial models"), # is not written in latex if float=F
            notes.append = T,
            notes = c("Standard errors are clustered at the level of countries."),
            covariate.labels = c("Concession (sum)",
                                 "Number of protesters (sum, log)",
                                 "Protester violence (sum)",
                                 "Demand: Policy change",
                                 "Demand: Institutional reforms",
                                 "Demand: Regime change",
                                 "State violence (sum)",
                                 "Physical violence index",
                                 "Democracy index",
                                 "Coup threat",
                                 "GDP pc (lag and log)",
                                 "Population (lag and log)",
                                 "Days since last protest",
                                 "Days since last protest (squared)",
                                 "Protest (lag)"),
            omit= c('year', 'gid', 'country'),
            omit.stat = c("logrank", "wald", "ll", "theta"),
            add.lines = list(c("Country dummies", "No", "Yes","Yes", "Yes"),
                             c("Year dummies", "No", "Yes","Yes", "Yes"),
                             c("Grid-year interaction", "No", "No", "No", "Yes"),
                             c("Number of countries", "18", "18", "18", "18")),
            label = "table:main-negbin-national",
            out = paste0("output/tables/Appendix_Table_varyingSpace",threshold,".tex"))
  

  
}


# Table B8 to B11: cluster SE at national level ---------------------------------------------------

threshold_list <- c(8,15,22,30)
for (threshold in threshold_list) {
  threshold1 = threshold -1
  
  load(paste0("data/navco_week_replication",threshold,".Rda"))
  

  nb1 <- glm.nb(protests_future ~ con_sum,
                data=navco_week)
  
  nb2 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  
  nb3 <- glm.nb(protests_future ~ con_sum + 
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  v2x_clphy + v2x_polyarchy +couprisk + nlights_calib_mean + gdp_log + pop_log +
                  #time
                  diff + diff2 + protests_past + pr_nr_c + 
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb4 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  #time
                  diff + diff2 + protests_past + pr_nr_c + 
                  factor(gid_year),
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))  
  
  
  
  # Table
  clustered_se <- list(
    sqrt(diag(vcovCL(nb1, vcov = vcovCL(nb1, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb2, vcov = vcovCL(nb2, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb3, vcov = vcovCL(nb3, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb4, vcov = vcovCL(nb4, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T)))
    
  )
  
  
  stargazer(nb1, nb2,nb3,nb4,
            se = clustered_se,
            type = "latex",
            float = F,
            float.env = "table",
            no.space = T,
            dep.var.caption =paste0("Sum of protests in the next ", threshold1, " days"),
            dep.var.labels = c(""),
            model.numbers = F,
            column.labels = c("1", "2", "3", "4"),
            title = paste0("Concessions and protests in the next ", threshold1, " days. Negative binomial models"), # is not written in latex if float=F
            notes.append = T,
            notes = c("Standard errors are clustered at the country level."),
            covariate.labels = c("Concession (sum)",
                                 "Number of protesters/protest (log)",
                                 "Protester violence (share)",
                                 "Demand: Policy change",
                                 "Demand: Institutional reforms",
                                 "Demand: Regime change",
                                 "State violence (share)",
                                 "Physical violence index",
                                 "Democracy index",
                                 "Coup threat",
                                 "Nightlights",
                                 "GDP pc (lag and log)",
                                 "Population (lag and log)",
                                 "Days since last protest",
                                 "Days since last protest (squared)",
                                 "Nr. of protest (lag)",
                                 "Protests in the rest of the country"),
            omit= c('year', 'gid', 'country'),
            omit.stat = c("logrank", "wald", "ll", "theta"),
            add.lines = list(c("Grid dummies", "No", "Yes","Yes", "Yes"),
                             c("Year dummies", "No", "Yes","Yes", "Yes"),
                             c("Grid-year interaction", "No", "No", "No", "Yes"),
                             c("Number of grids", "169", "169", "169", "169")),
            label = "table:main-negbin-national",
            out = paste0("output/tables/Appendix_Table_se_nat",threshold,".tex"))
  
  
}

# Tables B12 to B15: Share of concession ---------------------------------------------------

threshold_list <- c(8,15,22,30)
for (threshold in threshold_list) {
  threshold1 = threshold -1
  
  load(paste0("data/navco_week_replication",threshold,".Rda"))
  
  nb1 <- glm.nb(protests_future ~ con_share,
                data=navco_week)
  
  nb2 <- glm.nb(protests_future ~ con_share +
                  # control: protest characteristic
                  size_log+
                  protesterviolence_share + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence_share + 
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb3 <- glm.nb(protests_future ~ con_share +
                  # control: protest characteristic
                  size_log +
                  protesterviolence_share + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence_share + 
                  v2x_clphy + v2x_polyarchy +couprisk+ nlights_calib_mean + gdp_log + pop_log +
                  #time
                  diff + diff2 +protests_past + pr_nr_c + 
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb4 <- glm.nb(protests_future ~ con_share +
                  # control: protest characteristic
                  size_log +
                  protesterviolence_share + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence_share + 
                  #time
                  diff + diff2  +protests_past + pr_nr_c +
                  factor(gid_year),
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))  
  
  
  # Table
  clustered_se <- list(
    sqrt(diag(vcovCL(nb1, vcov = vcovCL(nb1, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb2, vcov = vcovCL(nb2, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb3, vcov = vcovCL(nb3, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb4, vcov = vcovCL(nb4, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T)))
  )
  
  stargazer(nb1, nb2,nb3,nb4,
            se = clustered_se,
            type = "latex",
            float = F,
            float.env = "table",
            no.space = T,
            dep.var.caption =paste0("Sum of protests in the next ", threshold1, " days"),
            dep.var.labels = c(""),
            model.numbers = F,
            column.labels = c("1", "2", "3", "4"),
            title = paste0("Concessions and protests in the next ", threshold1, " days. Negative binomial models"), # is not written in latex if float=F
            notes.append = T,
            notes = c("Standard errors are clustered at the level of grid cells."),
            covariate.labels = c("Concession (share)",
                                 "Number of protesters/protest (log)",
                                 "Protester violence (share)",
                                 "Demand: Policy change",
                                 "Demand: Institutional reforms",
                                 "Demand: Regime change",
                                 "State violence (share)",
                                 "Physical violence index",
                                 "Democracy index",
                                 "Coup threat",
                                 "Nightlights",
                                 "GDP pc (lag and log)",
                                 "Population (lag and log)",
                                 "Days since last protest",
                                 "Days since last protest (squared)",
                                 "Nr. of protest (lag)",
                                 "Protests in the rest of the country"),
            omit= c('year', 'gid', 'country'),
            omit.stat = c("logrank", "wald", "ll", "theta"),
            add.lines = list(c("Grid dummies", "No", "Yes","Yes", "Yes"),
                             c("Year dummies", "No", "Yes","Yes", "Yes"),
                             c("Grid-year interaction", "No", "No", "No", "Yes"),
                             c("Number of grids", "169", "169", "169", "169")),
            label = "table:main-negbin-national",
            out = paste0("output/tables/Appendix_Table_conShare",threshold,".tex"))
  
}



# Tables B16 to B19: concession dummy --------------------------------------------------------
threshold_list <- c(8,15,22,30)
for (threshold in threshold_list) {
  threshold1 = threshold -1
  load(paste0("data/navco_week_replication",threshold,".Rda"))
  
  nb1 <- glm.nb(protests_future ~ con_dummy,
                data=navco_week)
  
  nb2 <- glm.nb(protests_future ~con_dummy + 
                  # control: protest characteristic
                  size_log +
                  protesterviolence_dummy +  goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence_dummy +
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb3 <- glm.nb(protests_future ~ con_dummy +
                  # control: protest characteristic
                  size_log +
                  protesterviolence_dummy +  goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence_dummy +
                  v2x_clphy + v2x_polyarchy +couprisk + gdp_log + pop_log +
                  #time
                  diff + diff2 +protests_past + pr_nr_c +
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb4 <- glm.nb(protests_future ~ con_dummy +
                  # control: protest characteristic
                  size_log +
                  protesterviolence_dummy +  goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence_dummy +
                  #time
                  diff +  diff2 + protests_past + pr_nr_c +
                  factor(gid_year),
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  # Table
  clustered_se <- list(
    sqrt(diag(vcovCL(nb1, vcov = vcovCL(nb1, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb2, vcov = vcovCL(nb2, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb3, vcov = vcovCL(nb3, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb4, vcov = vcovCL(nb4, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T)))
  )
  
  
  stargazer(nb1, nb2,nb3,nb4,
            se = clustered_se,
            type = "latex",
            float = F,
            float.env = "table",
            no.space = T,
            dep.var.caption =paste0("Sum of protests in the next ", threshold1, " days"),
            dep.var.labels = c(""),
            model.numbers = F,
            column.labels = c("1", "2", "3", "4"),
            title = paste0("Concessions and protests in the next ", threshold1, " days. Negative binomial models"), # is not written in latex if float=F
            notes.append = T,
            notes = c("Standard errors are clustered at the level of grid cells."),
            covariate.labels = c("Concession (dummy)",
                                 "Number of protesters (log)",
                                 "Protester violence (dummy)",
                                 "Demand: Policy change",
                                 "Demand: Institutional reforms",
                                 "Demand: Regime change",
                                 "State violence (dummy)",
                                 "Physical violence index",
                                 "Democracy index",
                                 "Coup threat",
                                 "GDP pc (lag and log)",
                                 "Population (lag and log)",
                                 "Days since last protest",
                                 "Days since last protest (squared)",
                                 "Protest (lag)",
                                 "Protests in the rest of the country"),
            omit= c('year', 'gid', 'country'),
            omit.stat = c("logrank", "wald", "ll", "theta"),
            add.lines = list(c("Grid dummies", "No", "Yes","Yes", "Yes"),
                             c("Year dummies", "No", "Yes","Yes", "Yes"),
                             c("Grid-year interaction", "No", "No", "No", "Yes"),
                             c("Number of grids", "169", "169", "169", "169")),
            label = "table:main-negbin-national",
            out = paste0("output/tables/Appendix_Table_conDummy",threshold,".tex"))
  
}


# Tables B20 to B23: categorical concession --------------------------------------------------

threshold_list <- c(8,15,22,30)
for (threshold in threshold_list) {
  threshold1 = threshold -1
  load(paste0("data/navco_week_replication",threshold,".Rda"))
  
  nb1 <- glm.nb(protests_future ~ factor(con_cat),
                data=navco_week)
  
  nb2 <- glm.nb(protests_future ~ factor(con_cat) + 
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb3 <- glm.nb(protests_future ~ factor(con_cat) +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence +
                  v2x_clphy + v2x_polyarchy  +couprisk+ gdp_log + pop_log +
                  #time
                  diff + diff2 + protests_past +pr_nr_c +
                  gid_f + year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb4 <- glm.nb(protests_future ~ factor(con_cat) +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence +
                  #time
                  diff + diff2 + protests_past +pr_nr_c +
                  factor(gid_year),
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  # Table
  clustered_se <- list(
    sqrt(diag(vcovCL(nb1, vcov = vcovCL(nb1, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb2, vcov = vcovCL(nb2, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb3, vcov = vcovCL(nb3, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb4, vcov = vcovCL(nb4, cluster = navco_week[, "country_f"]),type = "HC1", fix = T, cadjust = T)))
  )
  
  
  stargazer(nb1, nb2,nb3,nb4,
            se = clustered_se,
            type = "latex",
            float = F,
            float.env = "table",
            no.space = T,
            dep.var.caption =paste0("Sum of protests in the next ", threshold1, " days"),
            dep.var.labels = c(""),
            model.numbers = F,
            column.labels = c("1", "2", "3", "4"),
            title = paste0("Concessions and protests in the next ", threshold1, " days. Negative binomial models"), # is not written in latex if float=F
            notes.append = T,
            notes = c("Standard errors are clustered at the level of grid cells."),
            covariate.labels = c("Neutral concession",
                                 "Non-material concession",
                                 "Material concession",
                                 "Full concession",
                                 "Number of protesters (log)",
                                 "Protester violence",
                                 "Demand: Policy change",
                                 "Demand: Institutional reforms",
                                 "Demand: Regime change",
                                 "State violence",
                                 "Physical violence index",
                                 "Democracy index",
                                 "Coup threat",
                                 "GDP pc (lag and log)",
                                 "Population (lag and log)",
                                 "Days since last protest",
                                 "Days since last protest (squared)",
                                 "Protest (lag)",
                                 "Protests in the rest of the country"),
            omit= c('year', 'gid', 'country'),
            omit.stat = c("logrank", "wald", "ll", "theta"),
            add.lines = list(c("Grid dummies", "No", "Yes","Yes", "Yes"),
                             c("Year dummies", "No", "Yes","Yes", "Yes"),
                             c("Grid-year interaction", "No", "No", "No", "Yes"),
                             c("Number of grids", "169", "169", "169", "169")),
            label = "table:main-negbin-national",
            out = paste0("output/tables/Appendix_Table_conCat",threshold,".tex"))
  
}



# Tables B24 to B27: zeroinflated negative binomial models ------------------------------------------------------------

threshold_list <- c(8, 15,22,30)
for (threshold in threshold_list) {
  threshold1 = threshold -1
  load(paste0("data/navco_week_replication",threshold,".Rda"))

  
  z1 <- zeroinfl(protests_future ~ con_sum,
                 dist = "negbin",
                 data=navco_week)
  
  z2 <- zeroinfl(protests_future ~ con_sum +
                   # control: protest characteristic
                   size_log +
                   protesterviolence + goal_policy + goal_reform + goal_regime +
                   # control: government characteristic
                   state_violence,
                 data=navco_week,
                 dist = "negbin",
                 control = zeroinfl.control(maxit = 1000, trace = 3))
  
  z3 <- zeroinfl(protests_future ~ con_sum +
                   # control: protest characteristic
                   size_log +
                   protesterviolence + goal_policy + goal_reform + goal_regime +
                   # control: government characteristic
                   state_violence + 
                   v2x_clphy + v2x_polyarchy +couprisk + nlights_calib_mean + gdp_log + pop_log +
                   #time
                   diff + diff2 + protests_past + pr_nr_c,
                 data=navco_week,
                 dist = "negbin",
                 control = zeroinfl.control(maxit = 1000, trace = 3))
  

  stargazer(z1, z2,z3,
            type = "latex",
            float = F,
            float.env = "table",
            no.space = T,
            dep.var.caption =paste0("Sum of protests in the next ", threshold1, " days"),
            dep.var.labels = c(""),
            model.numbers = F,
            column.labels = c("1", "2", "3", "4"),
            title = paste0("Concessions and protests in the next ", threshold1, " days. Zero-inflated negative binomial models"), # is not written in latex if float=F
            notes.append = T,
            covariate.labels = c("Concession (sum)",
                                 "Number of protesters (sum, log)",
                                 "Protester violence (sum)",
                                 "Demand: Policy change",
                                 "Demand: Institutional reforms",
                                 "Demand: Regime change",
                                 "State violence (sum)",
                                 "Physical violence index",
                                 "Democracy index",
                                 "Couprisk",
                                 "Nightlights",
                                 "GDP pc (lag and log)",
                                 "Population (lag and log)",
                                 "Days since last protest",
                                 "Days since last protest (squared)",
                                 "Protest (lag)",
                                 "Protests in the rest of the country"),
            omit= c('year', 'gid', 'country'),
            omit.stat = c("logrank", "wald", "ll", "theta"),
            add.lines = list(c("Grid dummies", "No", "No","No"),
                             c("Year dummies", "No", "No","No"),
                             c("Grid-year interaction", "No", "No", "No"),
                             c("Number of grids", "169", "169", "169")),
            label = "table:main-negbin-national",
            out = paste0("output/tables/Appendix_Table_zeroinf",threshold,".tex"))
}

# Tables B28 to B31: Poisson models -----------------------------------------------------------------

threshold_list <- c(8,15,22,30)
for (threshold in threshold_list) {
  threshold1 = threshold -1
  load(paste0("data/navco_week_replication",threshold,".Rda"))

  nb1 <- glm(protests_future ~ con_sum,
             family = poisson,
             data=navco_week)
  
  nb2 <- glm(protests_future ~ con_sum +
               # control: protest characteristic
               size_log +
               protesterviolence + goal_policy + goal_reform + goal_regime +
               # control: government characteristic
               state_violence + 
               gid_f + year_f,
             family = poisson,
             data=navco_week,
             control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb3 <- glm(protests_future ~ con_sum + 
               # control: protest characteristic
               size_log +
               protesterviolence + goal_policy + goal_reform + goal_regime +
               # control: government characteristic
               state_violence + 
               v2x_clphy + v2x_polyarchy +couprisk+ nlights_calib_mean + gdp_log + pop_log +
               #time
               diff + diff2  + protests_past + pr_nr_c +
               gid_f + year_f,
             data=navco_week,
             family = poisson,
             control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  nb4 <- glm(protests_future ~ con_sum +
               # control: protest characteristic
               size_log +
               protesterviolence + goal_policy + goal_reform + goal_regime +
               # control: government characteristic
               state_violence + 
               #time
               diff + diff2  + protests_past + pr_nr_c +
               factor(gid_year),
             data=navco_week,
             family = poisson,
             control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))  
  
  
  # Table
  clustered_se <- list(
    sqrt(diag(vcovCL(nb1, vcov = vcovCL(nb1, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb2, vcov = vcovCL(nb2, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb3, vcov = vcovCL(nb3, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb4, vcov = vcovCL(nb4, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T)))
    
  )
  
  
  stargazer(nb1, nb2,nb3,nb4,
            se = clustered_se,
            type = "latex",
            float = F,
            float.env = "table",
            no.space = T,
            dep.var.caption =paste0("Sum of protests in the next ", threshold1, " days"),
            dep.var.labels = c(""),
            model.numbers = F,
            column.labels = c("1", "2", "3", "4"),
            title = paste0("Concessions and protests in the next ", threshold1, " days. Poisson models"), # is not written in latex if float=F
            notes.append = T,
            notes = c("Standard errors are clustered at the level of grid cells."),
            covariate.labels = c("Concession (sum)",
                                 "Number of protesters (sum, log)",
                                 "Protester violence (sum)",
                                 "Demand: Policy change",
                                 "Demand: Institutional reforms",
                                 "Demand: Regime change",
                                 "State violence (sum)",
                                 "Physical violence index",
                                 "Democracy index",
                                 "Coup threat",
                                 "Nightlights",
                                 "GDP pc (lag and log)",
                                 "Population (lag and log)",
                                 "Days since last protest",
                                 "Days since last protest (squared)",
                                 "Protest (lag)",
                                 "Protests in the rest of the country"),
            omit= c('year', 'gid', 'country'),
            omit.stat = c("logrank", "wald", "ll", "theta"),
            add.lines = list(c("Grid dummies", "No", "Yes","Yes", "Yes"),
                             c("Year dummies", "No", "Yes","Yes", "Yes"),
                             c("Grid-year interaction", "No", "No", "No", "Yes"),
                             c("Number of grids", "169", "169", "169", "169")),
            label = "table:main-negbin-national",
            out = paste0("output/tables/Appendix_Table_poisson",threshold,".tex"))
  
}



# Tables B32 to B35: per regime type -----------------------------------------------------

threshold_list <- c(8,15,22,30)
for (threshold in threshold_list) {
  threshold1 = threshold%/%7
  load(paste0("data/navco_week_replication",threshold,".Rda"))

  nb1 <- glm.nb(protests_future ~ con_sum,
                data=navco_week)
  
  nb2 <- glm.nb(protests_future ~ con_sum +
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  
                  gid_f + 
                  year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  
  nb3 <- glm.nb(protests_future ~ con_sum*regime_reign_f + 
                  # control: protest characteristic
                  size_log +
                  protesterviolence + goal_policy + goal_reform + goal_regime +
                  # control: government characteristic
                  state_violence + 
                  v2x_clphy+v2x_polyarchy +couprisk + nlights_calib_mean + gdp_log + pop_log +
                  #time
                  diff +diff2  + protests_past + pr_nr_c + 
                  gid_f +  year_f,
                data=navco_week,
                control = glm.control(maxit = 1000, trace = 3, epsilon = 0.0001))
  
  
  
  # Table
  clustered_se <- list(
    sqrt(diag(vcovCL(nb1, vcov = vcovCL(nb1, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb2, vcov = vcovCL(nb2, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T))),
    sqrt(diag(vcovCL(nb3, vcov = vcovCL(nb3, cluster = navco_week[, "gid_f"]),type = "HC1", fix = T, cadjust = T)))
  )
  
  
  stargazer(nb1, nb2,nb3,
            se = clustered_se,
            type = "latex",
            float = F,
            float.env = "table",
            no.space = T,
            dep.var.caption =paste0("Sum of protests in the next week"),
            dep.var.labels = c(""),
            model.numbers = F,
            column.labels = c("1", "2", "3"),
            title = paste0("Concessions and protests in the next ", threshold1, " days. Negative binomial models"), # is not written in latex if float=F
            notes.append = T,
            notes = c("Standard errors are clustered at the level of grid cells."),
            covariate.labels = c("Concession (sum)",
                                 "Military regime",
                                 "Monarchy",
                                 "Party regime",
                                 "Personal dictatorship",
                                 "Number of protesters (sum, log)",
                                 "Protester violence (sum)",
                                 "Demand: Policy change",
                                 "Demand: Institutional reforms",
                                 "Demand: Regime change",
                                 "State violence (sum)",
                                 "Physical violence index",
                                 "Democracy index",
                                 "Coup threat",
                                 "Nightlights",
                                 "GDP pc (lag and log)",
                                 "Population (lag and log)",
                                 "Days since last protest",
                                 "Days since last protest (squared)",
                                 "Protest (lag)",
                                 "Protests in the rest of the country",
                                 "Concession (sum)*military regime",
                                 "Concession (sum)*monarchy",
                                 "Concession (sum)*party regime",
                                 "Concession (sum)*personal dictatorship"),
            omit= c('year', 'gid', 'country'),
            omit.stat = c("logrank", "wald", "ll", "theta"),
            add.lines = list(c("Grid dummies", "No", "Yes","Yes"),
                             c("Year dummies", "No", "Yes","Yes"),
                             c("Grid-year interaction", "No", "No", "No"),
                             c("Number of grids", "169", "169", "169")),
            label = "table:main-negbin-regime",
            out = paste0("output/tables/Appendix_Table_regime",threshold,".tex"))
  
}


# Figure B2 ---------------------------------------------------------------
threshold = 8
load(paste0("data/navco_week_replication",threshold,".Rda"))

# create new variable denoting the unique demands of protests per week in a protest wave
navco_week$goal <- vapply(strsplit(navco_week$camp_goals, ","), function(x) paste(unique(x), collapse = ","), character(1L))

navco_demand <- navco_week %>% 
  group_by(move_id) %>% 
  mutate(demand_change = ifelse(goal == lag(goal), 0, 1),
         demand_change = tidyr::replace_na(demand_change, 0),
         demand_nr = str_count(camp_goals, ",")+1,
         demand_nr_dummy = ifelse(demand_nr == 1, 0, 1)) %>% 
  ungroup() 

supp.labs <- c("No concession", "Concession")
names(supp.labs) <- c("0", "1")
navco_demand %>% 
  ggplot(aes(demand_nr_dummy)) +
  geom_histogram() +
  labs(
    x = "0: one demand; 1: over one demand",
    y = "n") +
  facet_wrap(~con_dummy, labeller = labeller(con_dummy = supp.labs, scales = "free")) +
  theme_lh
ggsave("output/figures/FigureB2.pdf", scale = 1, dpi = "retina", width = 16, height = 12, units = "cm")



# Figure B3 ---------------------------------------------------------------

threshold = 8
load(paste0("data/navco_week_replication",threshold,".Rda"))

# create new variable denoting the unique demands of protests per week in a protest wave
navco_week$goal <- vapply(strsplit(navco_week$camp_goals, ","), function(x) paste(unique(x), collapse = ","), character(1L))

navco_demand <- navco_week %>% 
  group_by(move_id) %>% 
  mutate(demand_change = ifelse(goal == lag(goal), 0, 1),
         demand_change = tidyr::replace_na(demand_change, 0),
         demand_nr = str_count(camp_goals, ",")+1,
         demand_nr_dummy = ifelse(demand_nr == 1, 0, 1)) %>% 
  ungroup() 

supp.labs <- c("No concession", "Concession")
names(supp.labs) <- c("0", "1")
navco_demand %>% 
  ggplot(aes(demand_change)) +
  geom_histogram() +
  labs(
    x = "Change in demands",
    y = "n") +
  facet_wrap(~con_dummy, labeller = labeller(con_dummy = supp.labs), scales = "free") +
  theme_lh

ggsave("output/figures/FigureB3.pdf", scale = 1, dpi = "retina", width = 16, height = 12, units = "cm")



